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We study via Monte Carlo simulation the dynamics of the Nagel-Schreckenberg model on a finite 
system of length L with open boundary conditions and parallel updates. We find numerically that 
in both the high and low density regimes the autocorrelation function of the system density behaves 
like 1— |t|/r with a finite support [— t, r]. This is in contrast to the usual exponential decay typical of 
equilibrium systems. Furthermore, our results suggest that in fact t = L/c, and in the special case of 
maximum velocity v max = 1 (corresponding to the totally asymmetric simple exclusion process) we 
can identify the exact dependence of c on the input, output and hopping rates. We also emphasize 
that the parameter r corresponds to the integrated autocorrelation time, which plays a fundamental 
role in quantifying the statistical errors in Monte Carlo simulations of these models. 
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1. INTRODUCTION 

The totally asymmetric simple exclusion process 
(TASEP) Q is a simple transport model, of fundamental 
importance in nonequilibrium statistical mechanics. In 
addition to its mathematical richness, it has applications 
ranging from molecular biology to freeway traffic. 

A TASEP consists of a chain of length L, with each 
site being either occupied by a particle or not, on which 
particles hop from left to right. See Fig.Q] If site i = 1 is 
vacant a particle will enter the system with probability 
a. If site i = L is occupied the particle will leave the 
system with probability (3. In the bulk of the system, a 
particle on site i will hop to site i + 1 with probably 1—p 
provided i + 1 is vacant, otherwise it remains at site i. 

TASEPs exhibit boundary-induced phase transitions, 
governed by the parameters a, (3 and p. In general, for a 
given p, there exist three possible phases, depending on 
a and /3: a low-density phase, a high-density phase, and 
a maximum-current (or maximum-flow) phase. 
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FIG. 1: A TASEP with L = 10 
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In the context of traffic models, it is most appro- 
priate to update all sites in parallel at each time-step. 
The stationary distribution of the TASEP with fully- 
parallel updates 0, Q is known exactly. (For reviews 
of the stationary properties of TASEPs with random se- 
quential updates see [3, The Nagel-Schreckenberg 
(NaSch) model Q is an important generalization of the 
parallel-update TASEP, in which particles can move up 
to Umax € N sites per time step. The NaSch model is 
generally considered to be the minimal model for traffic 
on freeways 0] ■ While many results are known rigorously 
for the TASEP, our understanding of the NaSch model 
and its further generalizations typically rely on numeri- 
cal simulation. This is particularly true of traffic network 
models, in which the NaSch model is often a component 
(see for example [clflcil]). 

In the current article we focus on dynamic 
(auto)correlation functions. The autocorrelations of the 
TASEP with random sequential update have been stud- 
ied in [111 Q3 an d display a separation of time scales 
between relaxation of local density fluctuations and col- 
lective domain wall motion. In particular, it was re- 
cently observed [l2| that the TASEP with random se- 
quential update exhibits non-trivial oscillations in the 
power spectrum of the system density, in the low and 
high density phases. In this article, we further elucidate 
the nature of these non-trivial oscillations, and demon- 
strate that they extend to the NaSch model generally. 
We emphasize that all the simulations performed in this 
work used fully-parallel updates, including our simula- 
tions of TASEP (which we view as the special case of the 
NaSch model with w max = 1). 
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1.1. Density autocorrelations 

The system density, n, which is simply the fraction 
of sites which are occupied, is an important quantity in 
many applications, including traffic modeling. The rela- 
tionship between density and flow is known as the fun- 
damental diagram in the traffic engineering literature. 
While the stationary-state expectation (n) of n is well 
understood for the NaSch model, and in fact known rig- 
orously for the TASEP, the dynamic behavior of n t is 
non-trivial. In this article we numerically study the au- 
tocorrelation function p n {t) '■= {(nont) — (n) 2 )/var(n) of 
the general NaSch model, and find a very simple form for 
its finite-size scaling. Up to very small corrections, our 
simulations show that in both the high and low density 
phases we simply have 



Pn(t) 



W 



1*1 < T, 
\t\ > T, 



(1) 



for some constant ral. 

The linear decay in ([T]) is in sharp contrast to the usual 
exponential decay typical of equilibrium systems. In fact, 
as discussed in section l2~2l there are good theoretical rea- 
sons to believe that p n (t) must ultimately decay exponen- 
tially on sufficiently long time scales, rather than exhibit 
the strictly finite support suggested by (JT|). However, as 
demonstrated by the simulations in sections [3] and 21 any 
corrections to the finite-support behavior displayed in ([!]) 
are extremely weak, and in practice (JT|) provides a very 
accurate approximation to the behavior of p n (t) through- 
out the low and high density phases. In particular, ([T]) 
provides a very good approximation to p n {t) for values 
of p relevant for traffic modeling. 

The Fourier series of p n (t) gives the power spectrum of 
n, and we note that taking the Fourier series of (H} does 
indeed produce oscillations as reported in [f2j]. Indeed, 
we have 
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The discussion in [12j focused on the case v max = 1 , with 
random sequential updates. However, our simulations 
show that ((T|) , and hence , hold more generally for the 
NaSch model with arbitrary v max . 

The specific form of the autocorrelation function 
has some interesting consequences for the design of Monte 
Carlo simulations. In particular, as discussed in sec- 
tion^ assuming the validity of ([1} we immediately have 
t = 2 Tint, n where r m t. n is the integrated autocorrelation 
time of n. The integrated autocorrelation time can be 
interpreted loosely as the number of time steps between 
"effectively independent" samples. It is therefore reason- 
able to conjecture that the parameter r should equal the 
amount of time it takes a fluctuation of the stationary 
state to traverse the system. If we let v denote the speed 
of such a fluctuation then we might reasonably expect 



that t = L/v. In section [3J we present numerical results 
that strongly suggest that in fact, for TASEP, we have 



r = L/\v c (a,P,p)\ 



(3) 



where v c (a, f3,p), the collective velocity [HUH, is known 
exactly. The results ([1} and ([3]) are consistent with the 
suggestions in [73 that the physical origins of the ob- 
served oscillations in the power spectrum of n are related 
to the time needed for a fluctuation to traverse the entire 
system. 

Furthermore, while no exact expression for 
v c (a, /3,p, Vmax) is known for the general NaSch model, 
the simulations presented in section [4] demonstrate that 
the scaling form ([3]) extends to general v max . In addi- 
tion, in the deterministic limit (p = 0) simple physical 
arguments produce an exact relationship between v c and 
u max which is in excellent agreement with the numerical 
results. 

The remainder of this article is organized as follows. 
In section [2 we briefly review some pertinent general 
theory relating to autocorrelations and then discuss some 
general consequences of ([T]). In section [3J we present our 
numerical evidence supporting ([T]) and {3J for TASEP, 
and also describe the exact expression for v c (a,p,p) in 
this case. We also explain relationship between ([TJ and 
© and the results presented in [f2j . In section |H we 
briefly review the definition of the NaSch model before 
presenting our numerical results for p n {t) in this case. 
Finally, we conclude in section [5] with a discussion. 

2. AUTOCORRELATIONS 

We begin by briefly recalling some standard definitions 
and results. Consider a Monte Carlo simulation of an 
ergodic Markov chain, and assume that sufficient time 
has passed that the system has reached stationarity. If 
one now measures an observable X at each time step 
one obtains a stationary time series X\ , X2 , ■ ■ ■ whose 
autocovariance function is defined to be 



C x {t) := (X X t )~(X o y< 



(4) 



The expectation (•) here is with respect to the stationary 
distribution, and we note that Cx(0) = var(A). The 
corresponding autocorrelation function is then defined as 



Px(t) 



C x (t) 
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(5) 



Finally, assuming Cx(t) to be absolutely summable, its 
Fourier transform defines the spectral density 
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C x (t), 



(6) 



The spectral density is closely related to the Fourier 
transform of the time series. Specifically, given any sta- 
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tionary time series X\ , X2 , • ■ • , Xt we can define its dis- 
crete Fourier transform to be 



t e 



(7) 



with oj = 2irm/T and m = 0, 1, . . . , T — 1. It is then 
straightforward to show [13] that for large T we have 



A!^)|V/xM+0(i 



(8) 



2.1.2. Exponential Autocorrelation time 

Typically, we expect that px{t) ~ exp(— i/"r cxp ) as 
t — > 00, which defines the exponential autocorrelation 
time T oxp . More precisely [TJj], one defines the exponen- 
tial autocorrelation time of observable X to be 



r eX p,x := hmsup — , 

Itl^oo log px{t) 



(14) 



and then the exponential autocorrelation time of the sys- 
tem as 



2.1. Autocorrelation times 



supr GXPi x, 
x 



(15) 



We now discuss the implications of the general form 
([T]) on two key time scales, the integrated autocorrelation 
time and the exponential autocorrelation time. 



2.1.1. Integrated Autocorrelation time 

From px(t) the integrated autocorrelation time is de- 
fined [151 as 
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00 

E 

t— — 00 



Px{t). 



If X denotes the sample mean of X\ , X2 , 
variance of X satisfies [13 



(9) 



. Xt then the 



var(A) ~ 2r inti x 



var(A) 
T : 



T 



(10) 



It is (|T0|) that accounts for the key role played by the 
integrated autocorrelation time in the statistical analy- 
sis of Markov-chain Monte Carlo time series. If instead 
of a correlated time series, one considers a sequence of 
independent random variables, then the variance of the 
sample mean is simply var(A)/X\ It is in this sense that 
Tint,x determines how many time steps we need to wait 
between two "effectively independent" samples. 

It can now be seen immediately from (jH]) that, as noted 
in the introduction, (JXJ) and ([3]) imply 



where the supremum is taken over all observables X. 
The autocorrelation time r oxp measures the decay rate 
of the slowest mode of the system, and it therefore sets 
the scale for the number of initial time steps to discard 
from a simulation, in order to avoid bias from initial non- 
stationarity. All observables that are not orthogonal to 
this slowest mode satisfy r cxp .x = 7 oxp . 

For the TASEP in continuous time, r oxp was computed 
analytically in [ll| [I?} using the exact Bethe Ansatz so- 
lution. In particular, it was found that r exp is O(l) with 
respect to L in the high and low density phases. We 
would expect the same behavior to hold generally for the 
NaSch model. 

However, if p n {t) were to have strictly finite support 
as claimed in ([1]), then we would have — \t\/ log p n (t) = 
for all \t\ > t, implying that r cxp! „ 7^ r cxp . This would 
then mean that n is orthogonal to the slowest relaxation 
mode, which seems implausible. We thus conclude that 
although ([T]) provides a very good approximation, p n {t) 
cannot actually have a strictly finite support. 



2.2. Finite-size scaling of p„(t) 

To obtain a more precise ansatz for p n {t) we therefore 
fix some ieN satisfying k < [t\ and set 



Pn{t) 
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Equation (| 13[) provides a very simple exact expression for 
7"int,n in terms of the physical parameters of the model. It 
is quite rare to have such an expression for a non-trivial 
model. 



Since we know empirically that ([T]) is a very good approx- 
imation, it must be the case that k/r ~ 1 as r — > 00. Let 
us then write r = k + e, where the only assumption we 
make regarding e is that e/r — > as t — > 00. Since 
the continuum limit of p(xt) should define a continuous 
function of x £ M we choose the parameter B by de- 
manding that 1 — \t\/r — _Be~' t ' Toxp when \t\ = k, which 
yields 



Pn(t) 



-(|i|-*0/r e: 



'/r, 



\t\ < k, 
\t\ > k. 



(17) 
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It is worth noting that the two expressions (Q} and (fl7|) 
lead to the same leading-order expression ([13"]) for r; nt „. 
Indeed, inserting (JTTj) into (j9]) we obtain 



2r, 



int.n 



e(1-e) + 



2t 



;l/r oxp _ I 



1 



(18) 



Since (e 1 /^ _ j/j-i = (!) for Tcxp = the terms 

arising from the exponential decay of p n (t) are 0(e) in 
the low and high density phases. 



3. TASEP 

We begin this section by comparing the power spec- 
trum found in [l2j with the Fourier transform of (H} . We 
then present the exact result for the collective velocity for 
TASEP, before presenting the results of our simulations. 



3.1. Power spectrum 

Let N denote the number of occupied sites in a TASEP 
system, and let N(w) denote the discrete Fourier trans- 
form of a particular time series Ni, N%, . ■ . , Nt, as de- 
fined in ©. The quantity := T(\N{lu)\ 2 ) is what 
[l2j refer to as the power spectrum of N. They find that 
for the continuous-time TASEP in the low-density phase 



T 



2v A 
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-Du> 2 L/v 3 



Llj 



(19) 



where A, D and v are parameters, which |12j set empiri- 
cally to v « 0.4, D w 20 and A « 1/500. 

We now attempt to compare (JT9J) with the corre- 
sponding result derived from (fTJ). From ((5J) we see that 
1 ^ 1 2 

/at(w) ~ ( iV(u;) ) as T — > 00, hence we should compare 
(]T9"]) with /jv(w) = L 2 f n (uj), where is computed 

via (JTJ. Although our empirical observations of the be- 
havior ([1} were made in the discrete time case of fully- 
parallel updates, (fTJ) can be interpreted as a well defined 
continuous function on R. In fact, the fully-parallel up- 
date rule becomes equivalent to the random sequential 
update in the limit e — > of rescaled variables 1 — p = e, 
a = ae and ft = fte. Here, a and ft are the usual injection 
and extraction rates of the TASEP in continuous time. 

To compare with the continuous time result (|19p . we 
compute / n (u;) via the continuous-time Fourier trans- 
form, so that ([1]) and Q predict 



= — — var(n) L 



1 



cos 



Luj 



(20) 



Now, since uj = 27rm/T, for sufficiently large T we have 
exp(Duj 2 L/v 3 ) ss 1. This is exactly the regime used 
by [II in their Fig. 3 (L = 1000 or L = 32000 and 
T = 10 6 ). Therefore, in this regime we can identify (|19[) 
with ([2D]) if v = \v c \ and 



.4 
D 



var(n) L. 



(21) 



Some remarks are in order. Firstly, for the determin- 
istic (p = 0) parallel-update TASEP, the static variance 
var(n) can be computed analytically from the known re- 
sults for the two-point function Q. In the low density 
phase it is given by 



a 



(1 



, 1 

3"! 



+ 0(L- 2 ), 



(22) 



and for the high density region a is replaced by ft. We 
expect that var(n) = 0(1/ L) would remain true when 
p > 0, and indeed for w max > 1 as well. In general, 
therefore, we expect the prefactor in (|2"0"|) to be 0(1) in 
L. 

Finally, we note that [12J fit (fT9|) to their data with 
a very small value of A/D. This small value follows 
from the fact that the numerical simulations in [l2T | were 
performed along the mean field line of the TASEP with 
random sequential update, where, theoretically, var(n) is 
identically zero. It is surprising that 12j were still able 
to extract a meaningful signal on this line. 



3.2. Collective velocity 

The stationary distribution of the TASEP with fully- 
parallel updates 0, Q is known exactly. In particular, if 
a < ft, 1 — yjp such TASEPs reside in a low-density phase, 
while for ft < a, 1 — */p a high-density phase results, with 
a = ft < 1 — y/p defining a coexistence line of the two 
phases (corresponding to a first order phase transition). 
For a, ft > 1 — yfp by contrast, the system resides in a 
maximum-current phase, in which the density is precisely 
1/2. 

The collective velocity |13| is the drift of the center of 
mass of a momentary local fluctuation of the stationary 
state, and is related to the current (flow) J and bulk 
density pt via v c — dJ(pb)/dpb- An exact expression for 
v c (a,ft,p) is available [2| for the case of parallel-update 
TASEP. If we define, for convenience, the function 



9(x,p) = 



(i-p){{i-xY-p) 

(1 - x) 2 +p(2x - 1) ' 



then 



v c (a,ft,p) 



-g(P,p), 



low density phase, 
high density phase. 



(23) 



(24) 



The negativity of the collective velocity in the high- 
density phase is simply due to the fact that it is the prop- 
agation of holes from right to left, rather than of particles 
from left to right, that is important in this phase. 

Using these exact expressions for v c the expression (J3]) 
now becomes 



\L/g(a,p), 
\L/g(P,p), 



a < ft, 1 - y/p, 
ft < a, 1 - ^/p. 



(25) 
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We note that for p — we have \v c \ = 1 identically 
throughout the high and low density regimes so that we 
simply have r = L in this case. We also note that in 
the low-density (high-density) phase r is independent of 
(3 (a). 



3.3. Simulations 

We now turn our attention to our Monte Carlo simu- 
lations. We simulated the parallel-update TASEP at a 
variety of values of a, /3 and p corresponding to both the 
low and high density phases, for system sizes L = 10 3 , 
5 x 10 3 and 10 4 . Each simulation consisted of lQ A L/v c 
iterations, with the first 10 3 L/v c time-steps discarded 
to ensure negligible bias due to initial non-stationarity 
(initially the system was empty). Assuming the validity 
of ([3|), this implies we generated 1.8 x 10 4 Ti n t,n samples 
of the stationary distribution in each simulation. 

For each simulation, we measured n at each iteration, 
and from the resulting time series we estimated the au- 
tocorrelation function p n {t) using the standard estima- 
tors |l5| . Fig. [5] shows a finite-size scaling plot of p„ (t) 
assuming the ansatz given by (fT]) with r = L, in the p = 
case. The agreement is clearly very good, and the sharp- 
ness of the cusp at t = L suggests that any corrections 
to the finite-support ansatz ([1} are very small. 



eluding the fundamental diagram) have cusps at p 
which are smoothed out for p > 0. 
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FIG. 2: Color online. Finite-size scaling plot of p n (t) for p = 
parallel-update TASEP in the high-density and low-density 
phases, for L = 10 3 , 5 x 10 3 , 10 4 and a variety of a, /3. 
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FIG. 3: Color online. Finite-size scaling plot of p n (t) for 
p = 0.25 parallel-update TASEP in the high-density and low- 
density phases, for L = 10 3 , 5 x 10 3 , 10 4 and a variety of a, /3. 
The choices of a, f3 shown correspond to four distinct values 
of v c providing strong evidence for the conjecture (|25[) . 
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FIG. 4: Color online. Finite-size scaling plot of pn(t) for 
p — 0.5 parallel-update TASEP in the high-density and low- 
density phases, for L = 10 3 , 5 x 10 3 , 10 4 and a variety of a, /3. 
The choices of a, j3 shown correspond to four distinct values 
of v c providing strong evidence for the conjecture (|25[) . 



Figs. [3] and |4] show finite-size scaling plots of p n (t) for 
p = 0.25,0.5, assuming the ansatz given by (fT]) and (1231) . 
There is again excellent data collapse, however we note 
that there is some noticeable curvature near the edge of 
the support, so that the sharp cusp present in the p = 
case becomes smoothed out somewhat for p > 0. As 
discussed in section 12.21 this does not affect the use of 
(TT5| for setting Monte Carlo error bars, but it would be 
interesting from a theoretical perspective to better un- 
derstand how this curvature depends on the model pa- 
rameters p, a, /3 and L (as well as v max ; c.f. the discussion 
in section [5]). We remark that many other quantities (in- 



4. NAGEL-SCHRECKENBERG MODEL 

An important generalization of the TASEP is the 
Nagel-Schreckenberg model @, in which each particle 
(vehicle) can move up to u max G N sites per itera- 
tion. Although the precise form of the phase diagram 
depends on w max , the NaSch model exhibits, in gen- 
eral, the same three qualitatively distinct phases as the 
TASEP [2l|. We now briefly review the dynamical rules 
defining the NaSch model. Suppose at time t £ N a 
vehicle with speed Vt G {0, 1, . . . , w max } is located on 
site Xt, and has headway (number of empty sites to its 
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right) equal to h t . Then the maximum speed this vehicle 
can safely achieve at the next time step is taken to be 
Vsafe = min(t> t + 1, Umax, ht), which allows for unit accel- 
eration provided the speed limit is obeyed and crashes 
are avoided. Provided u sa f e > 0, a random deceleration 
is then applied so that with probability p the new speed 
is v t +i = Wsafc - 1, otherwise v t +i = u safc . Finally, in 
the bulk of the system, the vehicle hops Vt+i sites to its 
right, so that Xt+i = x t + v t +\. All vehicles in the bulk 
of the system are updated in this way in parallel. The 
bulk dynamics clearly reduces to parallel-update TASEP 
when v max = 1. 

It remains to consider the boundary dynamics. We 
again wish to apply open boundary conditions, however 
choosing an appropriate implementation of such bound- 
ary conditions for the NaSch model is actually surpris- 
ingly subtle, and has been an active topic of research 
over recent years (l8l - [23| . In particular, it was argued 
in pll j that in order to observe the maximum-current 
phase when u max > 1 one needs to implement the inflow 
of vehicles into the system in a rather careful manner. 

Since our interest in the present context is confined to 
the high and low density phases however, we have chosen 
to implement the boundary conditions in the following 
simple way. We augment the system, which has sites 1 < 
i < L, with two boundary sites; one at i = and another 
at i = L + 1. With probability a a vehicle with speed 
Umax is inserted on site 0, and we immediately compute 
Usafe for this vehicle. If u sa f c > we move the vehicle to 
site u sa fe otherwise we delete it. The output is performed 
similarly. With probability 1 — /3 we insert a vehicle on 
site L+l, which then acts as a blockage to vehicles exiting 
the system. If the rightmost vehicle in the system has 
x t > L — Umax we define its new speed to be u sa f c and 
attempt to move the vehicle to site Xt+i = Xt + u sa f e . If 
Xt+i > L the vehicle is removed from the system. When 
Umax = 1 the above prescription reduces to the boundary 
rules for the simple TASEP described in section [TJ 

4.1. Simulations 

We now describe our simulations of the NaSch model 
as defined above. To our knowledge, no rigorous results 
are known for u c when u ma x > E However, for the deter- 
ministic case (p = 0) we expect that 

J Umax, low density phase, 
u c = < (26) 
1—1, high density phase, 

for any Umax, since in the low-density phase the deter- 
ministic movement of vehicles from left to right should 
control the dynamics, while in the high-density phase we 
expect that it is the movement of holes (traveling with 
speed 1) from right to left which is important. More gen- 
erally, we expect the form (|24j) to remain valid, but with 
an unknown function g, that will in general depend on 

Umax- 



Fig. [5] presents a finite-size scaling plot of p n (t) ob- 
tained by simulating the NaSch model with u max = 3 
and p = 0, with system sizes L — 10 3 , 5 x 10 3 and 10 4 
and a variety of values of a, (3 corresponding to both the 
low and high density phases. The data collapse is excel- 
lent, providing strong evidence for the ansatz obtained 
from (p}, ([3]), and ([26|). As for the case of p = when 
Umax = 1 we note the sharpness of the cusp at t = L/\v c \, 
again suggesting that any corrections to the ansatz ([1]) 
are very small. Each simulation performed consisted of 
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FIG. 5: Color online. Finite-size scaling plot of p n (t) for 
p = NaSch with t) max = 3 in the high-density and low- 
density phases, for a variety of choices of a, f3 and L. The 
exact value of v c is unknown in this case but here we chosen 
v c according to (12611 . 

10 4 L/|u c | iterations (with u c given by (1261) ). with the first 
10 3 L/|u c | time-steps discarded. The above simulations 
were also performed for u ma x = 5 with identical results. 

Finally, we also considered the case of u ma x — 3 with 
p = 0.25. For Umax > 1 and p > we are not aware of 
any exact predictions for u c , however it seems reasonable 
to conjecture that u c is independent of j3 {a) in the low 
(high) density phase. We therefore simulated the NaSch 
model with u ma x = 3, p = 0.25 and a = 0.25 at four 
different values of (3 > a, which should then correspond 
to a single value of v c . By considering a single value of 
v c we can still use a finite-size scaling plot of p n (t) to 
test the conjectures (P) and (J3]). Fig. [S] provides strong 
evidence to support their validity at u ma x > 1 and p > 0. 
By varying the value of |u c | used to produce the scaling 
plot of p n (t) so that the support edge lay at \v c \t/L = 1 
we obtained v c ~ 2.65. We remark that, assuming the 
validity of ([T]) and ^ , this method can be used as a way 
to obtain approximate values of |u c | when u ma x > 1 and 
p > 0. 



5. DISCUSSION 

We have studied the NaSch model in the low and high 
density phases via Monte Carlo simulation, and found 
that to a very good approximation the autocorrelation 
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\v c \t /L 

FIG. 6: Color online. Finite-size scaling plot of Pn(t) for 
p — 0.25 NaSch with i> ma x = 3 in the high-density and low- 
density phases, for a variety of choices of a, f3 and L. The 
exact value of v c is unknown in this case but here we have set 
v c = 2.65. 



function for the system density behaves as 1 — \v c t\/L 
with a finite support [— L/\v c \, L/\v c \], where v c is the 
collective velocity. For the case of f max = 1 an exact 
theoretical result is known for v c for all p <E [0, 1]. When 
f max > 1 no rigorous results for v c are known, however we 
conjecture that the when p = we simply have v c = v ma x 
in the low-density phase and v c = — 1 in the high-density 
phase. This result agrees with the exact result in the 
special case of v max = 1 and with numerical simulations 
for v max = 3,5. It seems reasonable to expect that it is 
valid for all w max for the deterministic NaSch model. 
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